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1-^ ' Abstract 

Let Q be an open, simply connected, and bounded region in R'', d > 2, 
and assume its boundary dQ is smooth. Consider solving an elliptic partial 
differential equation Lu = / over Q with zero Dirichlet boundary values. 
The problem is converted to an equivalent elliptic problem over the unit 
ball B, and then a spectral method is given that uses a special polynomial 
basis. With sufficiently smooth problem parameters, the method is shown 
' to have very fast convergence. Numerical examples illustrate exponential 

, convergence. 

§ : 1 INTRODUCTION 

o 



Consider solving the elliptic partial differential equation 

i«(s)^- ^ +^(s)«(s) = /(s), sencR'^ (1) 

with the Dirichlet boundary condition 

u(s) EE 0, sedn (2) 

Assume d > 2. Assume fl is an open, simply-connected, and bounded region in 
W^, and assume that its boundary dfl is several times continuously differentiable. 
Similarly, assume the functions 7(s), /(s), aij{s) are several times continuously 
differentiable over f2. As usual, assume the matrix A{s) = [aij(s)] is symmetric, 
and also assume it satisfies the strong ellipticity condition, 

€^^(s)C > coC^e, sen, ^eR" (3) 
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with Co > 0. Also assume 7(3) > 0, s G f7. 

In ^we consider the special region ft ~ B, the open unit ball in R''. We 
define a Galerkin method for with a special finite-dimensional subspace 

of polynomials, and we give an error analysis that shows rapid convergence of 
the method. In 331 we discuss the use of a transformation from a general region 
fl to the unit ball J5, showing that the transformed equation is again elliptic 
over B. Implementation issues are discussed in Sj4]for problems in K.^ and M'^. 
We conclude in ^ with numerical examples in and R'^ . 

The methods of this paper generalize to the equation 

3=1 

which contains first order derivative terms, provided the operator L is strongly 
elliptic. To do so, use the results of Brenner and Scott §§2.6-2.8], combined 
with the methods of the present paper. We have chosen to restrict our work to 
the more standard symmetric problem ([T]). 

There is a rich literature on spectral methods for solving partial differential 
equations. From the more recent literature, we cite [3], [16| . and [17j . Their 
bibliographies contain references to earlier papers on spectral methods. Our ap- 
proach is somewhat different than the standard approaches, as we are converting 
the partial differential equation to an equivalent problem on the unit disk or unit 
ball, and in the process we are required to work with a more complicated equa- 
tion. Our approach is reminiscent of the use of conformal mappings for planar 
problems. Conformal mappings can be used with our approach when working 
on planar problems, although having a conformal mapping is not necessary. 



2 A spectral method on the unit ball 



The Dirichlet problem has the following variational reformulation: Find 

u G Hq (r2) such that 



ds 



(4) 



/(s)u(s)ds, \fveH^{n) 



We define a spectral Galerkin method in this section for the special region 
n = B. In ^ we discuss the transformation of ([T]) from a general fl to an 
equivalent equation over the unit ball -B, a transformation that retains the 
ellipticity of the problem. In the remainder of this section, we replace fl with 
B. 
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Introduce the bilinear form 



dx, v,w e {B) 

and the bounded Uncar functional 

i{v)^ f /(x)«(x)dx, veH^iB) 
Jb 

The variational problem ^ can now be written as follows: find u E Hq (B) for 
which 

Aiu,v)^£{v), yveH^{B) (6) 
It is straightforward to show A is bounded, 

\A{v,w)\ < CA ||w||i 
CA = max||A(x)||2 + hiloo 

with ||-||-^ the norm of Hq (fl) and ||A(x)||2 the matrix 2-norm of the matrix 
A (x). In addition, we assume 

Aiv,v)>cM\l v^Hl[B) (7) 

This follows generally from ([3]) and the size of the function 7(x) over B\ when 
7 = 0, Ce = cq. Under standard assumptions on A, including the strong elliptic- 
ity in ([7]), the Lax-Milgram Theorem implies the existence of a unique solution 
u to ([6]) with 

Il"lli<-Il^ll 

Ce 

Denote by n„ the space of polynomials in d variables that are of degree < n: 
p (z n„ if it has the form 

2|<n 

with i a multi-integer, i = (ii, . . . , z^;), and \i\ = ii + ■ ■ ■ + id- Let denote our 
approximation subspace, 

A'„ = {(i-||x||^)p(x) |pen„} (8) 

with ||x||2 — xf + ■ ■ ■ + x^. The subspaces n„ and Xn have dimension 




A (w, w) 



J2 



9w(x) dw{x) 
dxi dxi 



+ 7(x)?j(x)w(x) 
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Lemma 1 Let A denote the Laplacian operator in W^. Then 

A : A-,, n„ (9) 

onto 

For a short proof, see [1] . 

The Galerkin method for obtaining an approximate solution to ([6]) is as 
follows: find u„ G Xn for which 

A{Un,v)^i{v), yveXn (10) 

The Lax-Milgram Theorem (cf. [3, §8.3], O §2.7]) implies the existence of m„ 
for all n. For the error in this Galerkin method, Cea's Lemma (cf. [3j p. 365], 
[3 p. 62]) implies the convergence of u„ to u, and moreover, 

||u-'Un||i<— inf (11) 

Ce v&X„ 

It remains to bound the best approximation error on the right side of this 
inequality. 

Given an arbitrary u G Hq {B), define w = — Am. Then w G [B) and u 
satisfies the boundary value problem 

-Au{P) = w{P), Pe B 
u{P) = 0, P edB 

It follows that 

u{P)= [ G{P,Q)w{Q)dQ, PeB (12) 
Jb 

For and R'^, the Green's function is defined as follows. 



d = 2: G{P,Q) = :^\og ^' 



2t: ^|r(P)-Qp 
d-3: G(P,Q) = --^ ' ^ ^ ^ 



47ri|P-Q| |P||r(P)-Q| 



(13) 



for P Q, Q G B, P ^ B. T{P) denotes the inverse point for P with respect 
to the unit sphere S'^'^ C M'*, 

T(rx) = ix, < r < 1, x G S"^"^ 
r 

Differentiate (1121) to obtain 



Vu{P)= [ [VpGiP,Q)]w{Q)dQ, PeB (14) 
Note that S/pG^P, •) is absolutely integrable over B, for all P e B. 
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Let Wn G n„ be an approximation of w, say in the norm of either C (B) or 
L2(B),andlet 

qn{P) = f G(P, Q)wn{Q) dQ, PeB 

JB 

We can show qn G Xn- This foUows from Lemma [Hand noting that the mapping 
in ([T2|) is the inverse of ([9]). 
Then we have 

u{P) - qn{P) = [ G{P, Q) [w{P) - Wn{Q)] dQ, PeB 

JB 

V [u (P) - 9„(P)] = / [VpG(P, Q)] [w{Q) ~ WniQ)] dQ, PeB 

JB 

The integral operators on the right side are weakly singular compact integral 
operators on (5) to (5) jJH chap. 7, §3]. This imphes 

\\U - qn\\l < c\\w - Wn\\a (15) 

By letting w„ be the orthogonal projection of w into n„, the right side will go to 
zero since the polynomials are dense in (B). In turn, this implies convergence 
in the [B) norm for the right side in pT|) provided u e Hq [B). 
The result 

inf — t)||i^O as n — + oo, u e (B) 

can be extended to any u e Hq (B). It basically follows from the denseness of 
Hq [B) in Hq (B). Let u e Hq (B). We need to find a sequence of polynomials 
{qn} for which g„||i — > 0. We know Hq {B) is dense in Hi (B). Given any 
k > 0, choose Uk e Hq (B) with ||it — Wfc||i < 1/fc. Then choose a polynomial Wk 
for which we have the corresponding polynomial q^ satisfying |jwfc — (7/c|| i < 1/fc, 
based on (fT5|) . [Regarding the earher notation, qk need not be of degree < fc.] 
Then \\u-qk\\i< 2/fc. 

To obtain orders of convergence, use (|15p and results on best multivari- 
ate polynomial approximation over the unit disk. For example, use results of 
Ragozin [15, Thm 3.4] or Yuan Xu [2U]. From [TS] we have the following. 

Theorem 2 Assume u e C^^"^ (^W) for some fc > 0. Then there is a polynomial 
qn e Xn for which 

\\u -qn\L<D (fc, d) n-' (n-l hlL,fc+2 + ^ ("''+'^ V")) (16) 

In this, 

ii«iL,.+2- E ll^'HL 

|j|<A;+2 

Lj{g,S)^ sup |g(x)-g(y)| 

|x-y|<5 

\i\=k+2 
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3 Transformation of the elliptic equation 



Consider the differential operator 

^.„s,^-i:A(«..(=)^) 

2.7 — 1 ^ ^ 



s e r2 c M'', w e (p) (17) 



which satisfies the elhpticity condition ([3]) with cq > 0. The operator M. is said 
to be elhptic on i?^ (fi). We want to transform the operator Al to one acting 
on functions u S C^{B) with B the unit baU in M'*. 

Assume the existence of a twice-differentiable mapping 



(18) 



and let ^ ^(^-^ M ^ B. Let 

onto 



J(x) = (i?$) (x) 



9(/3,,(x) 



xe B c 



denote the Jacobian of the transformation. As usual we assume J(x) is nonsin- 
gular on B, and furthermore 



min |det J(x)| > 



(19) 



Similarly, let K{s) = {D'^/) (s) denote the Jacobian of over fl. By differenti- 
ating the components of the equation 



* ($(x)) = X 



we obtain 



X($(x)) === (x), xe5 

This general approach is reminiscent of the coordinate transformations in [111 
Chap. 2] in which the mapping function is used in generating a mesh on a region 

n. _ 

For V e C^m, let 



D r- Tad 



v{x.) V ($ (x)) , X e B C 



and conversely, 
Then 



v{s) = u (* (s)) 

dv dv dipi (x) 
dxi dsi dxi 
'dipi (x) 



s e f7 C 



dsd dxi 
difd (x) 



1 ; 
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with the gradient VgW a column vector evaluated at s = $ (x). More concisely, 
Vx« (x) = J (x)^ V,v (s) , s = $ (x) (20) 

Similarly, 

Vs«(s) =if(s)TVxv(x), x = ^'(s) (21) 

Theorem 3 Assume the transformation $ satisfies or il9\} . Then for 

s = $(x), 

(^«)(^) ^-H^^wk^ t 4r fdet(J(x))a.,(x)^) (22) 



det (J(x)) .-^^ 9a;i \ 9x 



3 

I(x) X($(x))A($(x))Jf ($(x))'^ (23) 
= P«j(x)]fj=i 
Proof. Let w G (H) . Then 

{Mv){s)w{s)ds= I (7Wu)($(x))w($(x))det(J(x)) dx (24) 
n Jb 

On the other hand, using integration by parts we have 

d 



Jo. Jn j j^;^ osj osi 

f ^ ,,ai;($(x))9u;($(x)) , ,,,,,, 

(25) 

Using (dip, 

E (<j> (x)) ^-(;^w)^^(|jx)) ^ j^^^^^ ^^^^j ^^^^ j^^^^^ ^^^^j 

= [V^wi^)fK{<P (x))A ($ (x)) (x))^ [Vx^y(x)] 
= [VxiC(x)]^I(x) [V,w(x)] 
Using this to continue 



[ {Mv){s)w{s)ds^ [ [Vxw(x)]^A(x) [Vxw(x)]det(J(x)) dx 
Jn Jb 

^ - X E ^ (det (^(x)) «^,.(x)^) ^(x) dx (26) 
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Comparing ^24\ ) and i26]). and noting that w £ C§° (fl) is arbitrary, we have 

9 ^^„w7/„^^^ .^N^'^W 



(A^^;)(<i>(x))det(J(x)) = - ^ — f det (J(x)) S,,,(x) 



. . T 1 X 7 

which proves I122\) . ■ 

With this transformation, we can solve the Dirichlet problem over a general 
region by transforming it to an equivalent problem over the unit ball B. We 
can apply the Galerkin method to ^ by means of the transformation ((22l) . We 
convert ([1]) to the equation 

-E ^ (det(J(x))S,,,(x)^) +det(J(x))7(<I>(x)Rx) ^^7) 

= det(J(x))/($(x)) 

This system is also strongly elliptic. 

Theorem 4 Assume A{s), s e fi, satisfies and without loss of generality, 
assume 

det J(x) > 0, X e S 

Recall A (x) as defined by I123\) . Then A (x) satisfies the strong ellipticity con- 
dition 



Co = Co A* = Co min Amin(x) 



wiift Ai„in(x) the smallest eigenvalue of {x))'^K{^ (x)) (which equals the 
reciproc 
Proof. 



reciprocal of the largest eigenvalue of J (x) J (x) ) 



> CO {K^^) {K^C) - CO II 

In addition, 



|2 



||if(<i>(x))^e||2> Vin(x)||e|l2>A*||e||^ 

A* = min Amin(x) 

xS-B 

with Ajnin (x) i/ie smallest eigenvalue of {x))'^K{^ (x)),' c/. p. 
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4 Implementation 



Consider the implementation of the Galerkin method of Sl2]for the eUiptic prob- 
lem ([6]) over the unit ball B. We are to find the function u„ e Xn satisfying 
pU|) . To do so, we begin by selecting an orthonormal basis for n„, denoting it 
by {ipi, ■ ■ ■ , ^Pn}, with N = Nn = dimn„. Choosing an orthonormal basis is an 
attempt to have the linear system in (fTU)) be better conditioned. Next, let 



1 



1, 



(28) 



to form a basis for Xn- 
We seek 



Then ([TU)) becomes 



N 



(29) 



fc=i 



aV'fc(x) 97/'£(x) 
9a;j dxi 



+ 7(x)i/'fc(x)?Af(x) 



(30) 



= / /(x)7/^ax) dx, ^= 1,...,7V 



We need to calculate the orthonormal polynomials and their first partial deriva- 
tives; and we also need to approximate the integrals in the linear system. For 
an introduction to the topic of multivariate orthogonal polynomials, see Dunkl 
and Xu [7] and Xu [12]. For multivariate quadrature over the unit ball in R"^', 
see Stroud [18]. 



4.1 The planar case 

The dimension of n„ is 

7V„ = i (n + 1) (n + 2) (31) 

For notation, we replace x with (x, y). How do we choose the orthonormal basis 
{ipe{x, y)}£^i for n„? Unlike the situation for the single variable case, there are 
many possible orthonormal bases over B — D, the unit disk in M^. We have 
chosen one that is particularly convenient for our computations. These are the 
"ridge polynomials" introduced by Logan and Shepp [12| for solving an image 
reconstruction problem. We summarize here the results needed for our work. 
Let 

v„ = {Pe n„ : (P, Q) = vg e n„_i} 

the polynomials of degree n that are orthogonal to all elements of n„_i. Then 
the dimension of V„ is n -f 1; moreover, 

n„ = Vo ® Vi ® • • • ® Vn (32) 
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It is standard to construct orthonormal bases of each Vn and to then combine 
them to form an orthonormal basis of n„ using the latter decomposition. As 
an orthonormal basis of Vn we use 

1 TT 

fn,k{x,y) ^ —j=Un{x COS {kh) + y sin (kh)) , {x,y) e D, h= — — (33) 
y/TT n+1 

for fc = 0, 1, . . . , n. The function U„ is the Chebyshev polynomial of the second 
kind of degree n: 

Unit) = ^"^('^+^)^ t = cos6i, -1<<<1, 71 = 0,1,... (34) 
smt' 

The family {(ySn.fcjJ^^Q is an orthonormal basis of V„. As a basis of n„, we order 
{fn,k} lexicographically based on the ordering in ([55)1 and ([5^ : 

We}eLl = Wofii Vl,Oi V2,0, ■ ■ ■ , ^n,0, ■ ■ ■ , 'Pn,n} 



Returning to ([28)) . we define 

ipn,k{x,y) = (1 - a:^ - y^) (pn,k{x,y) (35) 

To calculate the first order partial derivatives of ipn,k{x, y), we need U„{t). The 
values of J7„(t) and J7„(t) are evaluated using the standard triple recursion 
relations 

Un+l(t) = 2tUn{t) - Un-l{t) 

Ui+,{t) = 2Un{t) + 2tui{t) - Ui_,{t) 

For the numerical approximation of the integrals in (|30p , which are over B 
being the unit disk, we use the formula 

J g{x,y)dxdy ^j^JldU ^) ^i^n (36) 

•"^ 1=0 m=0 ^ y ' / y ' 

Here the numbers uji are the weights of the (q + l)-point Gauss-Legendre quadra- 
ture formula on [0, 1]. Note that 



.1 9 

/ Pix)dx = ^p{ri)uji, 

Jo ,_o 



for all single- variable polynomials p{x) with deg (p) < 2q + 1. The formula ([36]) 
uses the trapezoidal rule with 2q + 1 subdivisions for the integration over B in 
the azimuthal variable. This quadrature is exact for all polynomials g € Il2q- 
This formula is also the basis of the hyperinterpolation formula discussed in [9] . 
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4.2 The three-dimensional case 

In M"^ , the dimension of n„ is 

Nn^ (^"3 ^) =i(n+l)(n + 2)(n + 3) 
Here we choose orthonormal polynomials on the unit ball as described in [7], 

= c„,,||x||™-^^pf ''"-^^■+^)(2||x|p - (^] , (37) 

j = 0, ...,[m/2j, /?==0,l,---,2(TO-2j), m = 0, l,...,n 

Here Cm.j = 2i~^^~^ is a constant, andpj*'''" ^■'+2)^ j ^ pj^^ g^j.g -j-j^g normalized 
Jabobi polynomials which are orthonormal on [— 1 , 1] with respect to the inner 
product 

(V, W)^ J (1+ i)""2j + |^,(i)y^(i) dt, 

see for example [1], [8]. The functions Sp^m-2j are spherical harmonic functions, 
and they are given in spherical coordinates by 



cos(f</>)T^^(cos6'), P even 
sm{^(j))T^ (cos 6*), P odd 



The constant c^^^ is chosen in such a way that the functions are orthonormal 
on the unit sphere 5^ in M"^: 



/ S'/3,fc(x)S'^^(x)d5 



The functions are the associated Legendre polynomials, see [10], [13]. Ac- 
cording to ([28]) we define the basis for our space of trial functions by 

V'mj,/3(x) = (1 - ||x||2).^™j,^(x) 

and we can order the basis lexicographically. To calculate all of the above func- 
tions we can use recursive algorithms similar to the one used for the Chebyshev 
polynomials. These algorithms also allow the calculation of the derivatives of 
each of these functions, see [8], [21] 

For the numerical approximation of the integrals in (|30p we use a quadrature 
formula for the unit ball B 

/•l p27T PIT 

g{x)dx^ / / g{r,e,(p)r'^sm{(j))d(l)dddrKiQg[g] 
Jo Jo Jo 

Qqig] "''S { ' arccos(O) I (38) 
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Here g{r, 9, cj)) = g{x) is the representation of g in spherical coordinates. For the 
integration we use the trapezoidal rule, because the function is 27r— periodic 
in 9. For the r direction we use the transformation 




where the v'/^ and ^fe are the weights and the nodes of the Gauss quadrature 
with q nodes on [— 1 , 1] with respect to the inner product 

{v,w)^J {1 + t)^v{t)w{t) dt 

The weights and nodes also depend on q but we omit this index. For the (j) 
direction we use the transformation 

/•I 

sm{<j))v{(j)) del) — / v{arccos{(j))) d(f> 

where the ujj and are the nodes and weights for the Gauss-Legendre quadra- 
ture on [—1, 1]. For more information on this quadrature rule on the unit ball 
in R3, see [Tg. 

Finally we need the gradient in Cartesian coordinates to approximate the 
integral in ([SO)) , but the function (pm,j,p{x) in ([57]) is given in spherical coordi- 
nates. Here we simply use the chain rule, with x = (x, y, z), 

^-v{r, 9, = ^vir, 9, 0) cos{9) sin(0) - ^v{r, P '''''^^^ 



dx dr 89 r sin(0) 

/ /, ,xCos(0)cos(0) 
aep r 

and similarly for and 

^ ay oz 



5 Numerical example 

Our programs are written in Matlab and can be obtained from the authors. 
Our transformations have been so chosen that we can invert explicitly the map- 
ping to be able to better construct our test examples. This is not needed 
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Figure 1: Images of (|4T|) . with a = 0.5, for lines of constant radius and constant 
azimuth on the unit disk. 



when applying the method; but it simplified the construction of our test cases. 
The elliptic equation being solved is 

Lu(s) = -Au + 7(s)u(s) = /(s), seSlCR'^ (39) 

which corresponds to choosing A ^ I. Then we need to calculate 

I(x)=X(<i>(x))X($(x))^ 



5.1 The planar case 

For our variables, we replace x e i? with (a;, y), and we replace s G with (s, t). 
Define the mapping ^ : B ^ Q hy {s,t) = ^ {x, y), 



s = X — y + ax 
t — X + y 



(41) 



with < a < 1. It can be shown that $ is a 1-1 mapping from the unit disk B. 
In particular, the inverse mapping ^! : ^ B is given by 



1 r 
a 
1 



-1 + y/l + a{s + t) 
at - (-1 + ^1 + a (s + t)j 



(42) 
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Table 1: Maximum errors in Galerkin solution u. 



n 


AT" 


II II 


7 

cond 


n 


AT 


II II 


cond 


n 
Z 





4.4i_n/ — i 


o.4z 


1 A 

14 


1 on 


y.yohi — D 


1/11 o 


3 


10 


4.21S- 1 


4.99 


15 


136 


3.03^;- 6 


165.8 


4 


15 


1.70E- 1 


9.27 


16 


153 


8.31£;- 7 


192.8 


5 


21 


9.63S- 2 


13.6 


17 


171 


2.09^;- 7 


222.1 


6 


28 


4.73S- 2 


20.7 


18 


190 


5.21^;- 8 


253.8 


7 


36 


1.88S- 2 


28.5 


19 


210 


1.42£:- 8 


287.9 


8 


45 


7.24£:- 3 


39.0 


20 


231 


3.53S- 9 


324.4 


9 


55 


2.79S- 3 


50.5 


21 


253 


7.58£;- 10 


363.4 


10 


66 


9.58£:-4 


64.7 


22 


276 


1.46£;- 10 


404.9 


11 


78 


3.20S-4 


80.4 


23 


300 


3.36£;- 11 


448.9 


12 


91 


9.67£;- 5 


98.6 


24 


325 


7.16£;- 12 


495.4 


13 


105 


3.01£:- 5 


118.7 


25 


351 


1.44£;- 12 


544.4 



In Figure [H we give the images in 17 of the circles r=j/10, j = l,...,10 and 
the azimuthal lines 9 — jtt/IO, j = 1, . . . , 20. 

The following information is needed when implementing the transformation 
from —Alt + 7u = / on f2 to a new equation on B: 

det (J) = 2(1 + ax) 

K 



A = KK 



det (J) A 



2(1 + ax) V -1 l + 2ax 
rp 1/1 ax 



2(1 + ax)^ \ '^^ '2,a X + 2aa:: + 1 



1 / 1 



ax 



1 + aa; V 2a x + 2ax + 1 



The latter are the coefficients for the transformed elliptic operator over B, given 
in dm). 

We give numerical results for solving the equation 

- Au(s,i)+e"~*'u(s,t) = /(s,t), {s,t) <E ft (43) 

As a test case, we choose 

M (s,t) = (l - - y^) cos (tts) (44) 

with (x,y) replaced using (H^ . The solution is pictured in Figure [21 To find 
f{s,t), we use P5)) and (HI]). We use the domain parameter a = 0.5, with fl 
pictured in Figure [T] 



14 



Figure 2: The true solution ([44]) 



Numerical results are given in Table [TJ The integrations in were per- 
formed with ([55)) : and the integration parameter q ranged from 10 to 30. We 
give the condition numbers of the linear system (|30p as produced in Matlab. 
To calculate the error, we evaluate the the numerical solution and the error on 
the grid 



The results are shown graphically in Figure El The use of a semi-log scale 
demonstrates the exponential convergence of the method as the degree increases. 

To examine experimentally the behaviour of the condition numbers for the 
linear system ([30)) . we have graphed the condition numbers from Table [T] in Fig- 
ure[4l Note that we are graphing Nn vs. the condition number of the associated 
linear system. The graph seems to indicate that the condition number of the 
system (pO|) is directly proportional to the order of the system, with the order 
given in (pij) . 



$ (xi J , j/i J- ) = $ {ri cos Oj , n sin 0j ) 
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Figure 3: Errors from Table [T] 



5.2 The three-dimensional case 

Here we define the niapping $ : _B ^ by (s, t, u) = $(a;, y, z), 

s = X — y + ax^ 
t = X + y 
u = 2z + hz^ 



(45) 



< a, 6 < 1, which is an extension of the mapping defined in (|41|) . The inverse 
mapping : il ^ _B is given by 

1 



a 
1 

a 

1 r 



-1 + \/l + a {s + t) 



at 



(-1 + + a(s + 0) 



-1 + Vl + bu 



In Figure [5] we show the image of the surface of B under As in the planar 
case, we also need 



r>$(a;,y, z) :=: J{x,y,z) 



1 + 2ax -1 
1 1 

2 + 2&Z 
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Figure 4: Condition numbers from Table [T] 



det(J(a;, y, z)) = 4(1 + ax){l + hz) 



and 



det[J{x,y,z))A{x,y,z) 

= det( J (x, y, z)K{x, y, z)K'^{x, y, z) 

( 1 



4(1 + ax){\ + hz) 



2(l + aa;)2 2(l + aa;)2 
ax 1 + oa; + 2c?x^ 



2(l + aa;)2 2(l + aa;)2 




4(l + 6z)2 / 



Again, these are the coefBcients for the second order term for the transformed 
equation on i?, given in (j22p . We give numerical results for solving the equation 

— Aw(s, i, u) + e'*^*w(s, i, u) = /(s, i, u), (s, i, m) g 

and for our test case we choose 

v{s, t, u) = sin Q(s - t)^ ■ (1 - ||^'(s, t, u)f ) 

where the second term guarantees the Dirichlet boundary conditions on fl. 
Numerical results are given in Table [2] 
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Figure 5: Image of ((45|) from two different angles, with a=0.7, b=0.9, for lines 
of constant ip and 9 on the sphere. 



The integrations in (|30p were performed with ([38]); and the integration pa- 
rameter q was chosen as q = n + 2. Numerical experiments indicate that a 
larger q does not change the results significantly. The condition numbers for 
the system ([50]) were again calculated with Matlab. An estimation for the 
error in the maximum norm was calculated on the grid given by 




/ ^ sin ( A^) cos \ 

^ sin (^^) sin (f tt) 
V ^cos(i-^) / 



1,...,20, j = l, 



,40. 



The error for the Galerkin method is shown in Figure [S] and the development 
of the condition number is shown in Figure [T] Again the numerical experiment 
seems to indicate an exponential convergence of the method and a linear growth 
of the condition numbers with respect to the number of degrees of freedom Nn 
of the linear system ((30|) . 



Additional Remarks. We present and study a spectral method for the Neu- 
mann problem 

-Au(s) -|-7(s)u(s) = /(s), sef^CM'' 

in a forthcoming paper. We are also investigating the behaviour of the con- 
dition number for the linear system ([50)1 associated with our spectral method, 
attempting to prove that it has size 0(iV„), consistent with the numbers shown 
in Tables [Hand [2 
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Table 2: Maximum errors in Galerkin solution u. 



n 


AT 


M M 


7 

cond 


1 


4 


A AO 77" 1 

4.98-C/ — 1 


1.5 


2 


10 


1.99i/ — 1 


3.0 


Q 
O 


on 
ZD 


i. ( o£/ — 1 


0. / 


4 


35 


8.22£;-2 


11.0 


5 


56 


2.18i;-2 


17.1 


6 


84 


1.34£; - 2 


27.1 


7 


120 


5.95E - 3 


39.4 


8 


165 


i.eoi; - 3 


55.9 


9 


220 


4.85£; - 4 


75.8 


10 


286 


2.56£;-4 


100.2 


11 


364 


1.44£; - 4 


128.9 


12 


455 


7.85i; - 5 


162.4 


13 


560 


4.19£;-5 


200.6 


14 


680 


2.33i; - 5 


244.0 



Our earlier numerical examples use given 

onto 

chosen to be nontrivial and illustrative. In general, however, when given a 
smooth mapping 

ifi-.S^dil, 

onto 

it may not be clear as to how to extend tp to ^ over B. In some cases there is 
an obvious choice, as when O is an ellipsoid. We are investigating schemes to 
produce continuously differentiable extensions $ which satisfy 

min |det J(x)| > 
and for which J(x) is easily computable. 
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Figure 6: Errors from Table [2] 
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